Gravitational dynamics of an infinite shuffled lattice: early time evolution and 

universality of non-linear correlations 



T. Bacrtschiger 

Dipartimento di Fisica, Universita "La Sapienza" , P.le A. Mow 2, 1-00185 Rome, Italy, 
& ISC-CNR, Via dei Taurim 19, 1-00185 Rome, Italy. 



OO 
O 
O 

(N 
X 



o 

a 



c 

o 
o 



> 

On 
(N 



O 
O 



X 



M. Joyce 

Laboratoire de Physique Nucleaire et de Hautes Energies, UMR 7585, 
Universite Pierre et Marie Curie — Paris 6, 75252 Paris Cedex 05, France. 

F. Sylos Labini and B. Marcos 
"E. Fermi" Center, Via Panisperna 89 A, Compendio del Viminale, 1-00184 Rome, Italy, 
& ISC-CNR, Via dei Taurim 19, 1-00185 Rome, 



In two recent articles a detailed study has been presented of the out of equilibrium dynamics 
of an infinite system of self-gravitating points initially located on a randomly perturbed lattice. 
In this article we extend the treatment of the early time phase during which strong non-linear 
correlations first develop, prior to the onset of "self-similar" scaling in the two point correlation 
function. We establish more directly, using appropriate modifications of the numerical integration, 
that the development of these correlations can be well described by an approximation of the evolution 
in two phases: a first perturbative phase in which particles' displacements are small compared to the 
lattice spacing, and a subsequent phase in which particles interact only with their nearest neighbor. 
For the range of initial amplitudes considered we show that the first phase can be well approximated 
as a transformation of the perturbed lattice configuration into a Poisson distribution at the relevant 
scales. This appears to explain the "universality" of the spatial dependence of the asymptotic 
non-linear clustering observed from both shuffled lattice and Poisson initial conditions. 



PACS numbers: 05.40.-a, 95.30.Sf 



I. INTRODUCTION 



Structure formation in the universe is currently ad- 
dressed primarily using numerical simulations of purely 
self-gravitating particle systems, with initial configura- 
tions generated by displacing the particles slightly from 
a perfect lattice (see e.g. 0> §])• The physics of the 
strongly non- linear regime of the observed evolution is, in 
detail, very poorly understood. Progress in understand- 
ing would be very useful in providing analytical guidance 
for numerical simulations, and in particular better con- 
trol on their precision in representing the relevant con- 
tinuum limit. In a scries of recent papers |}], ^] we have 
studied a reduced version of the full cosmological prob- 
lem, considering a very simple class of randomly per- 
turbed lattices as initial conditions, and evolution in a 
static universe 1 . 

One of our primary results is that, despite the simpli- 
fications, the system we study has qualitative behavior 
very similar to that observed in the more complex cos- 
mological simulations (correlated perturbations in initial 
conditions, expanding universe). Notably the evolution is 
clearly "hierarchical" (i.e. structures build up at succes- 
sively larger scales driven by the linearized fluid theory 



growth of the initial perturbations), and asymptotically 
"self-similar" (i.e. the time dependence of the two point 
correlation function is given by a simple scaling of the 
spatial variable which can be inferred from the linearized 
fluid theory) 2 . The functional form of the spatial de- 
pendence of the non-linear correlation function is, on the 
other hand, just as in the cosmological simulations, a 
fundamental quantity characterizing the bresults which 
is determined numerically, but not currently understood 
(i.e. not predicted analytically or semi-analytically). We 
have noted in ||, Q], however, that this form emerges, to 
a good approximation, in our simulations prior to the 
asymptotic bscaling regime, in the preceding transient 
phase in which strong non-linear correlations first de- 
velop. In this paper we extend and detail our analysis 
of this phase. We show in greater detail that the emer- 
gence of the observed non-linear two body correlations 
can be very well approximated by modeling the evolu- 
tion as constituted of two subsequent phases, with an 
abrupt matching from one to the other. During the first 
phase the particles evolve as described by a perturbative 



1 See also || 
conditions. 



for earlier studies of evolution from these initial 



2 To avoid any possible confusion we note that both these terms are 
used here with meanings different to those commonly ascribed to 
them in condensed matter physics. In the latter context both are 
associated with invariance proper ties of the spatial correlations 
under spatial rescalings (see e.g. |l l| ) . Such properties are not 
implied by their use in the present context. 
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analytical approximation we have introduced and studied 
in 0, ||; in the second phase the particles evolve under a 
force coming solely from their nearest neighbor (NN). 

We relate our work here to previous work along similar 
lines concerning evolution from Poissonian initial condi- 
tions [|[ ||, D . In this case it has been shown explicitly || 
that the emergence of the first strong non-linear correla- 
tions can be very well accounted for by an approximation, 
at the relevant early times, in which the full gravitational 
force on each particle is truncated to that due only to 
its initial NN 3 . In the case of "shuffled lattice" (SL) ini- 
tial conditions, which we consider in this work, such an 
approximation is not generically good: when the typical 
displacement of a particle is small compared to the lattice 
spacing, the high degree of symmetry gives that the force 
on a typical particle is the sum of comparable contribu- 
tions from many particles. In this regime, however, we 
can describe the evolution very well by a simple pertur- 
bative approximation, which has been developed fully in 
[Q, ||. This latter approximation breaks down, roughly, 
when particles start to approach one another, which is 
precisely when one expects a NN approximation for the 
force may become appropriate. We show here that this 
is indeed the case, and that an abrupt switch between 
the two phases gives a very good approximation to the 
evolution. Further this model allows us to explain the 
fact that the observed form of the non-linear correlations 
in our simulations is independent of the amplitude of the 
initial shuffling, and the same as that observed from Pois- 
son initial conditions. This is the case because, for the 
range of amplitude of the initial perturbations we use, the 
evolution in the first phase brings the system to a distri- 
bution with correlation properties which, at the relevant 
scales, are essentially those of a Poisson distribution. 

The main interest of our results here is that they give 
a semi-analytical understanding of the origin of the form 
of the observed non-linear two point correlations for this 
class of initial conditions, which are qualitatively simi- 
lar to those used in cosmological type simulations. As 
remarked in |^| the form of this early time correla- 
tion function coincides, to a very good approximation, 
with that which is also observed in the asymptotic scal- 
ing regime attained by the system at longer times. This 
suggests strongly (but does not prove) that the physical 
mechanism leading to the former correlations, which we 
identify here, is also that which gives rise to the latter cor- 
relations. In the context of cosmological simulations such 
a conclusion, if appropriate also for that case, would be 
very important for the following reason. In this context 
the results derived from the numerical simulations which 
are physically relevant are those which are representa- 
tive of the Vlasov-Poisson (VP) limit of the simulated 



particle system. The mechanism we describe here for the 
generation of the non- linear correlations is, on the other 
hand, clearly not representative of this limit: the effect 
of interactions with single NN particles are precisely of 
the kind which are discarded in the VP limit, which is a 
mean-field approximation. Therefore, if the form of the 
non-linear correlations in the long time evolution turns 
out actually to be determined in such a phase, this form 
would not be representative, as required, of the VP limit. 
We discuss this point a little further in our conclusions, 
and suggest numerical tests which could be performed 
to determine whether the long-time behavior is indeed 
linked to the early time mechanism we study here. 

The paper is organized as follows. In Sec. || we discuss 
some of the relevant properties of the initial conditions 
of our simulations. In the next section we give the de- 
tails of the simulations considered here and summarize 
briefly the main relevant results of || In Sec. [II 



3 The importance of NN interactions at early times starting from 
Poisson initial conditions has also been discussed previously by 



Saslaw. See 



and references therein. 



present in detail the two-phase model which captures the 
essential elements of the formation of the first non-linear 
correlations. We give results here also of numerical sim- 
ulations. Finally in Sec. IV we discuss the results and 
draw our main conclusions. 



II. STATISTICAL PROPERTIES OF INITIAL 
CONDITIONS 

As in jj| we study evolution from initial conditions in 
which the particles are at rest and located at the sites of 
a perfect simple cubic lattice subjected to random uncor- 
rected displacements. We adopt the same notation, de- 
noting by p(u) the probability density function (PDF) for 
the displacements, and by u(R) the displacement of the 
particle originally at lattice site R. The variance of the 
PDF is denoted by A 2 , and the dimensionless variance 
5 2 = A 2 /£ 2 , where £ is the lattice spacing. The param- 
eter 6 we refer to as the normalized shuffling parameter. 
As discussed in Q| , for the case of purely gravitational 
interactions, the system is completely characterized, in 
the infinite volume limit, by the single parameter i5. The 
Poisson distribution corresponds to the limit in which 
each particle's position is completely randomized in the 
infinite volume, i.e., S = oo. 

The precise details of the different initial conditions 
of which the evolution is studied numerically below, are 
summarized in Tab. Q. The PDF p(u) used for gener- 
ating the displacements is constant in a cube of side 
2A around the origin (and with sides parallel to the 
axes of the lattice), and zero elsewhere. The (arbitrary) 
choice of units is as in ||] , giving that the dynamical time 
Td V n = l/y^JrUpo is equal (where po is the mass density). 

The first four simulations (SL64 to SL16) are the same 
ones analyzed in S. As explained there (see Sec. Ill 
of H) the values of 6 have been chosen so that, in our 
units of length, 5 2 l b is constant. This gives (see ||] for 
details) an amplitude of the power spectrum at small k 
(i.e. k <C l^ 1 for 5 < 1) which is equal in all Simula- 
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Name 


N i/a 


L 


I 


A 


5 


m/m§4 


SL64 


64 


1 


0.015625 


0.015625 


1 


1 


SL32 


32 


1 


0.03125 


0.0553 


0.177 


8 


SL24 


24 


1 


0.041667 


0.00359 


0.0861 


18.96 


SL16 


16 


1 


0.0625 


0.00195 


0.03125 


64 


SL64b 


64 


1 


0.015625 


0.0012 


0.0768 


1 


P64 


64 


1 


0.015625 


00 


oo 


1 



TABLE I: Details of the initial conditions studied in this pa- 
per, and the numerical parameters used in the simulations. 
N is the number of particles in the cubic box of side L, and 
m is the particle mass. 



tions. With time in units of Td yn this means that, in 
the long wavelength fluid limit, the systems are identical 
initially, and evolve identically. SL64b differs only from 
SL64 in the value of 8, i.e., they are two simulations with 
identical values of the parameters characterizing the fi- 
nite numerical representation of the infinite systems, but 
with different 8. 

To characterize the correlation properties of the distri- 
butions we will use the same quantities as in §, §: the 
reduced two point correlation function £(r), the power 
spectrum P(k) (which is related to £(r) by a Fourier 
transform 4 ), the NN PDF ui(r). We refer the reader to 
H for the precise definitions of these quantities. We will 
also consider the stochastic properties of the force, which 
we characterize using P(F), the PDF for the modulus of 
the force F. 

While 8 — oo corresponds exactly to the Poisson dis- 
tribution, one expects any SL with 8 > 1 to approximate 
the correlation properties of a Poisson distribution up to 
a scale of order A = 81. Indeed, if 8 > 1, the effect of the 
short distance exclusion of the underlying lattice should 
disappear and the particles are, to a good approxima- 
tion, randomly placed in a volume of order A 3 . This can 
be seen explicitly for the power spectrum, of which the 
exact analytical expression may be written ^ [llj in the 
form 



P(k)^- + \p(k)\ 2 A(k) 
n 



(1) 



where no is the mean particle density and p(k) is the 
characteristic function of p(u) (i.e its Fourier transform 
normalized so that p(0) = 1). The function A(k) depends 
only on the initial unperturbed lattice distribution (and 
not on the shuffling). For any simple form of the PDF 
for the shuffling (such as the top-hat one considered here 
and in ||, or, e.g., a Gaussian PDF as used in Q) p(k) 
decreases toward zero for fcA > 1, giving that the power 



4 The power spectrum is the Fourier transform of £(r), which dif- 
fers from what we refer to here as the "correlation function" by 



spectrum tends to the Poissonian value (given by 1/no). 
Thus for wave- numbers k larger than of order l/(5£) the 
power spectrum converges to this Poissonian behavior. 
We refer the reader notably to Fig. 2 of §3], which show 
the power spectrum for the four initial conditions SL16, 
SL24, SL32 and SL64. 



A. Nearest neighbor distribution 

We will consider below often the NN PDF, and it is 
useful to know its form in the initial conditions just de- 
scribed. For the Poissonian limit 5 = oo it is straightfor- 
ward to show analytically ITTI] that it is 



uipir) = Airrior exp 



(2) 



which gives an average distance between NN Aq ~ 

— 1/3 

0.55n /J = 0.551 Since the NN distribution character- 
izes the small scale properties we expect, following our 
discussion above, that this expression will be a good ap- 
proximation for 8 > 1. For 8 <C 1, on the other hand, 
one may show 5 that 



£w(r) 



(28)'- 



where 



8(28 + x) [2Sx 
[8 2 - x 2 + 8x] 
x [28x 2 - 
f(8-x 




±S 3 

3 U 



4S 2 x + fx 



(3) 



if a; € [-28, -8} 

if x e ] - 8, 0] 
if x e ]0, 8} 
otherwise, 



corresponding to average distance between NN Ao = £ — 
(86827/80640) A. 

In Fig. [I] we show the behavior of u)(r) for most of 
the SL studied here. For the SL with very small shuf- 
fling — SL16 or SL24 — this PDF is strongly peaked 
around the average distance between NN (which is ap- 
proximately equal to £), and in very close agreement with 
the analytical approximation given by Eq. (||) . For SL32 
(with 8 = 0.177) a small discrepancy with this approx- 
imation is visible, while for SL64 (with (5 = 1) the NN 
PDF is, as expected, in very good agreement with that 
in the Poisson case given by Eq. (0). 



a delta function singularity at r = 0. See 



or 



The derivation of the expression given is straighforward, but te- 
dious. For a given shuffling of a particle and its six NN, one 
must determine exhaustively the different combinations, and as- 
sociated probabilities, which lead to a given NN separation. The 
approximation 6 <C 1 is used in taking the inter-particle separa- 
tions to linear order in S. 
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FIG. 1: NN PDF of the SL considered in this paper (the name 
of each SL is indicated above the corresponding curve). For 
SL32, SL24 and SL16, the function (^) is shown for compari- 
son. In the insert panel we show enlarged the SL64 (SL) case 
together with the behavior of 0J p (r) for a Poisson distribution 
(P) with same number density [i.e. Eq.(^)]. 

B. Force distribution 

The PDF of the modulus of the force W(F) is a useful 
quantity in our analysis. Notably if the forces on particles 
arc dominated by that coming from their NN the simple 
relation 

W(F)dF = uo(r)dr (4) 

must hold. For the case of a Poisson distribution the 
analytical expression for W(F) was first given by Chan- 
drasekhar jL2| . It is proportional to the so-called Holtz- 
mark distribution (see |IT| ] for the explicit result and a 
simple derivation). In Fig. ^ we show a plot of the full 
(Holtzmark) distribution, and the PDF inferred if only 
NNs contribute Wnn(F)cIF — uj{r)dr. The domination 
of the NN in the force is clearly seen at stronger values 
of the force. The relation is not valid at weaker values of 
the force as these correspond to the (rare) particles for 
which the force picks up comparable (and possibly can- 
celing) contributions from more than one particle. Note 
that the tail of the PDF at large F decays in proportion 
to F~ 5 / 2 , which means that the variance (i.e. second 
moment) of the PDF is infinite. 

In Jl3| we have studied in detail the statistical proper- 
ties of the force in an SL. As one might expect, one can 
show that the force PDF is very well approximated by 
that of a Poisson distribution when the typical displace- 
ment is larger than the inter-particle spacing, i.e., 5 > 1. 
At small values of the displacements, on the other hand, 
the force PDF is very different to that in a Poisson distri- 
bution, decaying much more rapidly at large values of the 
force: the strong forces due to NN are completely absent 
as the typical particle feels a comparable effect from its 
six NN when the configuration is close to a perfect simple 
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FIG. 2: Holtzmark distribution and the PDF inferred if only 
NNs contribute, i.e., WNN(F)dF = cu(r)dr. The agreement is 
very good in the strong field limit where W(F) ~ F~ 5 ^ 2 . At 
weak field the PDF due to the NN has a sharp cut-off while 
the Holtzmark distribution shows a more gentle decrease (see 
discussion in |l|). 

cubic lattice. More precisely, for a top-hat PDF of the 
displacements (as used here), the functional behavior of 
the PDF at large F changes qualitatively, from exponen- 
tial decay to a F~ 5 / 2 power law decay, at 6 = 0.5. For 
6 > 0.5 the amplitude of the latter tail is lower than that 
in the Poisson, with this difference becoming negligible 
as 8 increases to of order unity. 



III. NUMERICAL SIMULATIONS 

In this section we first report results of numerical sim- 
ulations in which the initial conditions given in Tab. || 
are evolved under the mutual self-gravity of all particles. 
As such simulations have already been reported in detail 
in |j, |3] (see also |5|, |j| for Poisson initial conditions and 
some SL), we restrict ourselves to a very brief summary 
with an emphasis on the points which are relevant here. 
We then report results of a new set of simulations de- 
signed to validate our model of the early time evolution 
by a direct numerical integration of the appropriate two 
phase approximation. 



A. Full gravity 

As in [ pi El we have used the publicly available code 
Gadget [|14[ |l5|| to evolve the system under gravity, 
modified only by a small scale regularization of the po- 
tential below r = e. This softening parameter e is taken 
here in all simulations to be e = 0.00175L (i.e. in all sim- 
ulations significantly smaller than the initial average dis- 
tance between NN) . We have performed the same checks 
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r 




FIG. 3: "Collapse plot" of £(r,t) (absolute value) for the P64 
initial conditions: for each time t > 1 we have rescaled r so 
that £(r, t) = 1 at ro, where £(ro,i = 1)- The behavior of the 
numerical fit given in |3J is also shown for comparison. 



FIG. 4: Evolution of the function R s (t) in P64 (points) com- 
pared with its prediction from linear fluid theory, R s (t) oc 
exp[(2/3)t/rdyn]- 



as discussed in ||, || for the independence of our results 
to this choice. 

In H we have found that the SL initial conditions con- 
sidered here lead, from the time significant non-linear 
correlation first develops at small scales, to an evolution 
in which the correlation function can be approximated 
by 

£(r,t)«3(r/R.(t)) , (5) 

where R s (t) is a time dependent length scale, and a sim- 
ple functional fit to S(r) is given in ||. For sufficiently 
long times — after a transient time of which the duration 
increases as the value of S decreases — R s (t) follows very 
well the behavior predicted by a simple analysis based on 
the linearized equations for the system approximated as 
a pressure- less self-gravitating fluid. As described in 
for such a system with power spectrum of density fluctu- 
ations at small wave-numbers P(k) ~ k n one obtains: 

R a (t) <x exp C_H —\ . ( 6 ) 

V 3 + n T dyn J 

As reported in Q SL initial conditions indeed produce 
the predicted asymptotic time dependence, correspond- 
ing to n = 2 in this formula. In Fig. |^ we show the same 
analysis of the evolved Poisson initial conditions P64, us- 
ing the same numerical fit to 5(r) as found in || for the 
SL initial conditions. In Fig. ffl is shown the associated 
temporal evolution of R s (t), and a fit to the theoretical 
fluid behavior, given by Eq. (JsJ) with n = 0. The agree- 
ment, after a short transient, is as good as that observed 
in (| for the SL. 

The conclusion which follows is thus that, while the 
temporal behavior of the scaling depends on the value of 
5, the functional form of the spatial dependence in the 
non-linear correlation function appears to be the same for 



all initial conditions. This is what we refer to as "univer- 
sality" of the non- linear correlations in this context. 

We note that, just as underlined for the SL initial con- 
ditions in (|, m, it is also true for the Poisson initial 
conditions that the spatio-temporal scaling of the two 
point correlation function, given by Eq. (|5[), is a good 
approximation well before the asymptotic scaling behav- 
ior, given by Eq. ([j]), sets in. While the dynamical model 
we present below is valid only in this first phase (i.e. 
that prior to the asymptotic regime), it is thus natural 
to hypothesize that the form of the asymptotic correla- 
tion function is in fact determined in this phase. We will 
discuss this hypothesis further in our conclusions, and in 
particular, how it could be tested for. 

The second essential result about the development of 
non-linear correlations which we recall is the following. 
In all these simulations (both Poisson [|j and SL || [|] ) 
we observe that the relation 

oj(r) dr = (l- J w ( s ) ds ^j ■ 4irr 2 n (1 + f (r)) dr , (7) 

holds to a very good approximation, from the time that 
significant non-linear correlations first develop until a 
time of order a dynamical time later. As explained in 
H (see also JTH), it is valid if all but the two point cor- 
relations are trivial. It is thus natural to interpret its 
observed approximate validity for the correlations which 
develop in the first phase of non-linearity to indicate that 
these correlations develop predominantly as a result of 
the two body clustering of NN pairs of particles. For the 
case of the Poisson initial conditions it has been shown 
explicitly in that this interpretation is correct: by 
integrating from the initial conditions with only forces 
between initial NN pairs, the evolution of correlation is 
well described up to approximately one dynamical time 
when non-linear correlations have developed up to a scale 
of order I. 
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FIG. 5: Schematic representation of the evolution of two point 
correlations [specifically f(r) + 1] from SL initial conditions 
with small S. During the first phase (1) the initial anti- 
correlation (i.e. exclusion) at small scales is destroyed, as 
particles evolve in the regime described by PLT. In the sec- 
ond phase (2) large positive correlations are created at small 
scales, up to roughly the initial lattice spacing. The forces re- 
sponsible for these correlations are predominantly those of ex- 
erted by NN pairs on one another. In the subsequent evolution 
(3), when this approximation is no longer valid, the regime 
of positive correlations grows in a self-similar way, larger and 
larger scales becoming non-linear with time. 



B. Two phase model evolution 



FIG. 6: The two-point correlation function at the times 
tmax = 1, 2.5, 3.5, 4.5 rdyn for the different initial conditions 
as indicated, in both the full gravity simulations (thick lines) 
and the simulations of the two phase model described in the 
text (thin lines). The corresponding transition times t* are 
the optimal ones given in Eq. (^). For clarity the x axis 
has been rescaled for each initial condition (as otherwise the 
curves are, to a very good approximation, all superimposed, 

cf. IB). 



tion emerges, leads us to consider the approximation of 
the early time evolution in which one abruptly matches 
a PLT phase onto a NN dominated phase, i.e., 



For SL initial conditions, with small 5, the approxima- 
tion of forces as NN dominated is, as we have discussed 
above, not valid at early times. In this limit, however, 
we have developed in |Q ||] an analytical perturbative 
approach, which at linear order gives a very good ap- 
proximation to the dynamical evolution. The treatment 
involves simply a Taylor expansion of the force between 
particles in their relative displacements from their ini- 
tial lattice positions R, and thus breaks down when the 
latter become equal to the initial separation of the parti- 
cles. In H the precision of the linearized approximation 
has been explored in detail for the SL initial conditions 
(and others). For the evolution of the average relative 
distance between NN, the approximation turns out to 
be very good until this quantity becomes quite close to 
the initial lattice spacing, i.e., until when many particles 
come close to their NN. We refer to this approximation 
as particle linear theory (PLT) as it is simply a general- 
ization for particles of an analogous standard treatment 
for the self-gravitating fluid (in the Lagrangian formal- 
ism, leading to the so-called Zeldovich approximation). 
We do not detail further the implementation of this PLT 
approximation here as a succinct summary may be found 
in 0j , and a very complete discussion in || . 

The fact that PLT is observed to work very well up to 
close to the time of NN domination, and the observation 
that Eq.(^) is valid when significant non- linear correla- 



• Phase 1: From t = up to a time t» particles in 
the system evolves according to PLT. 

• Phase 2: For t > t* particles evolve only subject 
to the gravitational attraction of their NN at the 
time t*. 

While the approximation used in the first phase is good 
for the whole system, and in particular describes well 
the evolution of correlations at any distance, the second 
phase will only be valid approximately in describing cor- 
relation at some sufficiently small scale, and for suffi- 
ciently short times. A schematic representation of the 
evolution of correlations is given in Fig. |^. 

We have implemented the above two phase evolution 
numerically on the set of initial conditions given in Tab. |. 
We have taken the time at which we match the approx- 
imations, £*, as a free parameter and adjusted it to best 
fit the evolution of the correlations in the full gravity 
simulations in the phase when strong non-linear correla- 
tion first emerge 6 . We find that for each initial condition 
there is indeed a choice of the time which gives such 



6 In the second (NN) phase we use the same numerical value of 
e = 0.00175 as in the full gravity simulations (and the same 
functional form of the smoothing as in GADGET). 



a fit, to a very good approximation over a range of am- 
plitudes from £(r) ~ 10 2 down to considerably less than 
unity. Results are shown in Fig. for the optimal times 

U ~ 0, 0.5, 1.5 and 3.0 (8) 

for SL64, SL32, SL24 and SL16 respectively (with time 
in units of Tdyn)- The results here are given at the times 
t max , which are the approximate times at which we ob- 
serve the evolution under NN interactions to lead to cor- 
relations beginning to deviate from (i.e. lag behind) those 
in the full gravity simulations: 

t max = 1, 2.5, 3.5 and 4.5. (9) 

For times t > t max the modified simulations stop evolving 
significantly (as one is then simply seeing the averaged 
effect of the periodic motion of many NN pairs). In con- 
trast the full gravity simulations continue to evolve clus- 
tering from the collective motion of larger scales which 
has been completely neglected in the second phase of the 
approximation. 




FIG. 7: The NN PDF in the distributions evolved to the 
times tMN given in Eq. |lo| for the different SL initial condi- 
tions indicated. Also shown is the result for a Poisson distri- 
bution [Eq. (H)], and a line indicating its small scale behavior 
ujp(r) oc r 2 . Units on the x (y) axis have been divided (mul- 
tinliedl bv the numerical value of the lattice snacm? f. in each 



C. Transition to NN domination 

We now consider in more detail the essential time scale 
t*, which we determined numerically above. Given our 
discussion of and motivation for the two phase model, we 
might expect it to correspond to the time at which PLT 
breaks down and the forces on particles become typically 
dominated by that due to their NN. As we will now ex- 
plain it corresponds in fact to a time somewhat shorter 
than this. 

An approximate characterization of the time tNN at 
which NN forces dominate can be given numerically by 
studying the relation between the NN P DF uo {r) and the 
force PDF W{F). As explained in Sec. |n| NN forces 
dominate when the relation Eq. (^) holds, for large val- 
ues of the modulus of the force F. While in the initial 
condition SL64 it already holds to a very good approxi- 
mation, in the distributions obtained by evolution of the 
SL with smaller S we find that it becomes good at the 
times 



t NN « 2.0, 3.0 and 4.2 



(10) 



for SL32, SL24 and SL16 respectively. Fig. ^ shows the 
NN PDF in each case at these times, and Fig. || the 
PDF of the force modulus. The quantities have been nor- 
malized by the characteristic length/force scale in each 
case to give the PDF for the corresponding dimension- 
less quantity: for the NN PDF we have normalized the 
radial distance to the initial lattice spacing £ in each case, 
and for the force PDF we have normalized to the force 



(2ml)/- 



dyn 



Gm 2 /£ 2 (i.e. the force between two par- 



ticles at the characteristic distance £). We have plotted 
the two quantities separately to show the fact that at the 
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FIG. 8: PDF of the force at the times tNN given in Eq. |l^ 
for the different SL initial conditions as indicated. The force 
has been rescaled as described in the text. Also shown is the 
PDF for the P64 initial conditions, as well as the theoretical 
result for an infinite PDF (i.e. the Holtzmark distribution). 



times tNN we find in all cases good agreement at small 
separations/strong forces with the corresponding results 
for a Poisson distribution. This will be important in our 
discussion below of the origin of the "univers ality " of the 
observed correlations. As discussed in Sec. II B the re- 
lation Eq. 



is very accurately followed for a Poisson 
distribution (cf. Fig. ||), and so it follows implicitly from 
these figures that it is indeed a good approximation here 
in the regime of strong forces. The very clear differences 
in the force PDF at smaller values of the force are a re- 



8 



flection of the real difference in the fluctuations at larger 
scales: the larger amplitude at smaller values of the force, 
for a force PDF normalized as we have done, is a reflec- 
tion of the fact that the fluctuations become more sup- 
pressed at large scales as 6 decreases. In the NN PDF, 
on the other hand, such differences due to large scales 
are not present (because it characterizes the small scale 
properties). 

It is clear, from Eqs. (||) and (|To|), that t* does not 
correspond, even approximately, to tNN, but rather is 
shorter by more than a dynamical time. Likewise it is 
shorter than the estimated time of breakdown of PLT. 
Indeed we find that at tNN we have in all cases that the 
root square variance of the displacements of NN parti- 
cles, normalized in units of the appropriate £, is ~ 0.5 
which corresponds to the criterion for breakdown of PLT 
found in ||. At i*, on the other hand, we measure the 
values 

6 rs 0.18, 0.12, 0.11 (11) 

for SL32, SL24 and SL16 respectively 7 . 

This difference, and the approximate value of £», can 
in fact be understood quite easily as follows. An approx- 
imation to the evolution which is clearly better than our 
two phase approximation at all times is that in which we 
use the expansion of PLT to approximate all the forces on 
a particle except that due to the particle which becomes 
its NN in the second phase. Indeed we would expect such 
an approximation — let is call it PLT+NN — to be very 
good for the whole regime we are treating. If we now 
consider our two phase treatment (PLT for t < t*, NN 
for t > as an approximation to PLT+NN, it is not 
difficult to understand why the i* which makes the ap- 
proximation optimal is of order a dynamical time smaller 
than t^N- The reason is that the equation of motion for 
the displacements in PLT reduces, for time scales up to 
of order a dynamical time, to a simple ballistic approx- 
imation 8 . Thus when we turn on the NN interaction at 
time t — f* we follow, for of order a dynamical time, 
the PLT+NN evolution. If we reduce i* further we will 
deviate from PLT+NN more because PLT is poorly ap- 
proximated; if we increase i* we will loose precision by 
excluding the full NN contribution to the force. 



These values are, as expected, also in very good agreement with 
those predicted by PLT. A more appropriate characterisation of 
the breakdown of PLT is in fact given by considering the variance 
of the relative displacements. The difference with respect to the 
simple (one-point) variance is in fact negligible here for reasons 
that will be explained below. 

It is straightforward to show that one obtains u(R, t) = u(R, 0)+ 
v(R, 0) t, where v(R, 0) is the velocity of the particle associated 
with the lattice site R at the (arbitrary) initial time t = 0, by 
expanding to linear order the full PLT expression for the evolu- 
tion of the displacements (see H, or [pj) in powers of e n (k)t/r t j yn 
where e n (k) are numbers of order unity specifying the eigenval- 
ues of the dynamical matrix for gravity on the lattice jd H] . 



Note that PLT "comes for free" in this way for a time 
after t* only because we match the velocities at t*. A 
simple check on the above explanation is provided by 
doing a simulation in which we reset the velocities to 
zero at t*. We have done this and find indeed that the 
match to the full gravity evolution is considerably less 
good, and that, particularly for the SL with smaller 8 no 
choice of can produce a good fit. We will return in the 
next section to discuss further the role of the velocities. 



D. Origin of universality of non-linear correlations 

The simple two phase model thus describes quite well 
the emergence of the non-linear correlations at early 
times for the range of initial conditions studied. We now 
turn to the fact that these correlations are approximately 
the same in all cases, i.e., the two point correlation func- 
tions, as a function of radial separation, agree quite well 
(and indeed agree well with these quantities as found in 
| and |,1). 



1. Small scale correlation properties at tNN 

The explanation for this "universality" of the cluster- 
ing is clearly suggested by the results shown above in 
Figs. |?] and || in all cases the evolution gives at t « tNN 
a point distribution with correlation properties very sim- 
ilar to those of the Poisson distribution, at the scales 
relevant to the development of clustering in the follow- 
ing phase. Indeed the force PDF, at stronger values of 
the force corresponding to the particles which will clus- 
ter most rapidly in the subsequent time, follows very 
closely that in the Poisson distribution. Thus the clus- 
tering then develops as in the latter distribution giving 
the same correlation properties. These are simply those 
which emerge, as described in detail in j^] , when pairs of 
particles with initial separations given by the NN PDF 
top(r) in the Poisson distribution fall on one another. 

Why are the correlation properties so similar at these 
scales to those of a Poisson distribution at the time ijvw? 
We have emphasized in Sec. || that a Poisson distribution 
is approximated to increasingly large scales as the ampli- 
tude S of the random uncorrelated shuffling increases. It 
follows that if the evolution described by PLT is to a 
good approximation simply an amplification of the ini- 
tial displacements, with only very weak correlation, the 
transformation to a "universal" Poisson initial condition 
at the relevant scales at t « tNN results. To quantify 
whether this is indeed the case we can consider, e.g., 

fEr (u(0) • u(Rjvjv)) 

CNN{t) = r^T 

where Rnn ar e the lattice vectors of the six particles of 
the simple cubic lattice closest to the origin (where there 




is a particle). The ensemble average is the average over 
realizations of the random initial conditions of the SL 9 . 

In Fig. H is shown cjvjv(i) a s calculated exactly in PLT. 
We see that, for the time scales over which we use PLT (at 
most about four dynamical times), the correlation which 
develops between displacements at the relevant (small) 
scales is indeed weak (cat at < 20%). 

2. The limit 6^0 

Fig. [I shows that if, instead, our initial conditions had 
S sufficiently small so that t^N were greater than a few 
dynamical times, the approximation of weak correlation 
of the displacements at small scales would become pro- 
gressively worse as 6 decreases. As a result the basis 
for the approximate universality in the subsequent evo- 
lution would also be expected to become a progressively 
poorer approximation. The evolution of the displace- 
ments in PLT is simply a sum over the appropriately 
evolved eigenmodes of the displacements fields in the 
corresponding linear approximation to the inter-particle 
force. The behavior we observe here at long times is 
a result, as discussed in detail in ^, [§], of the fact 
that in this regime the small spread in the eigenvalues of 
the modes of the displacement field becomes important. 
The modes with slightly faster growth become arbitrar- 
ily dominant, leading to the very specific correlation of 
displacements described by these modes. For arbitrarily 
long times, i.e., for arbitrarily small initial S, one there- 
fore obtains a distribution with correlation properties at 



When this average is performed all six approximate NN are 
equivalent so that it is in fact sufficient to evaluate cjviv(i) for a 
single neighbor [and drop the sum and factor of 1 /6 in Eq. (02)] . 



all scales very different to that which could form from the 
Poisson distribution 10 . We thus conclude that the "uni- 
versality" we observe in our numerical simulations is a 
good approximation in the range of SL initial conditions 
with small, but not very small 11 , <5. 



3. Role of velocities attistN 

In the above discussion we have neglected the role of 
the velocities: in SL32, SL24 and SL16 non-zero veloc- 
ities have developed at tNN which make the full initial 
conditions at this time different from those in SL64 and 
P64 (with vanishing velocities at t^N = 0). The PDF 
W v (v) of the modulus of these velocities as measured in 
the different simulations at this time are shown in Fig. [Tc| . 
We have normalized for convenience in units of 

Gm , „. 

^Vt 1 (13) 

which is the velocity gained by a particle initially at rest 
when it reduces its separation by one half to a particle 
initially at distance t. 

We observe that the different PDF of the velocities 
agree quite well, which means that the full (space and 
velocity) initial conditions at tNN are indeed very sim- 
ilar in these simulations. Compared to SL64 and P64, 



Specifically, as S — > the evolution will always be dominated at 
the time when PLT breaks down bv the most rapidly growing 
eigenmode. In an SL lattice (see JH, H]) this eigenmode is one in 
which adjacent infinite parallel planes fall towards one another. 
Such small initial 5 are very difficult to simulate numerically 
because of the precision required. 
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FIG. 11: Behavior of s(t) [as defined in Eq. (|14|)] in SL32, 
SL24 and SL16. 



however, the difference in velocities is a priori signifi- 
cant: their magnitude is not small, but of order unity, 
in the units chosen, which are characteristic for the next 
stage of free fall of NN particles which leads to the cor- 
relations. Given that at this time tjyN we expect the 
velocity of particles to be, on average, oriented towards 
their NN (since t^jy is significantly larger than i*), and 
that in the approximated Poisson distribution the aver- 
age NN distance is 0.55^, the distribution at t^N should 
be well approximated as one in which pairs of NN fall on 
one another, but starting at an earlier time. 

This is further illustrated and quantified by Fig. [ll], 
which shows the temporal evolution of the quantity 



e(t) = 



(v ■ i-jvjy) 

|v||rArjv| 



(14) 



for t > 0, where v is the velocity of a given particle (at 
time t > 0), and rjvjv is the vector pointing towards its 
NN (at the same time t), and the average is over all the 
particles. The evolution of this quantity is qualitatively 
very similar in all three simulations: after a slow rise from 
an initial non-zero value a rapid decrease sets in at a time 
close to the estimated tjvjv in each case [cf. Eq. (10)]. The 
characteristic time for this, roughly exponential, decay 
of the correlations corresponds well to At = £ max — t nn 
[where £ max is as estimated in Eq. (j^)]. We thus see, 
as anticipated, that there is significant correlation of the 
direction of velocity with the NN direction at time tjviv- 
While At « 1 for SL64 and PL64, we have At « 0.5 
for the three simulations shown here, with a very similar 
behaviour of the correlation function s(t) in this phase. 
Thus SL32, SL24 and SL16 all lead to similar non-linear 
correlations as those in SL64 and P64, but in a shorter 
time due to this correlation of velocities with the NN 
direction acquired before i^Ar. 

The behavior of s(t) observed here can be understood 
in greater detail in the model we have described. The ini- 



tial non-zero, approximately constant, value is a result of 
the fact that at sufficiently small times PLT can be well 
approximated by its fluid limit. In this case 0, ^ the 
displacement of each particle off its lattice site is simply 
amplified in time, so that the function s(t) is indepen- 
dent of time and equal to the expression in Eq. ( |l4| ) with 
v replaced by u, the initial displacement of the particle 
from its lattice site (giving s(0) ~ 0.6). The slow increase 
of correlation is due to the difference between PLT and 
this fluid limit. The decrease from about £jvjv signals 
that pairs of NN particles have now begun to cross one 
another, giving a contribution with the opposite sign to 
s(t). At sufficiently long times a given particle's motion 
is finally no longer oriented with the direction of its NN, 
as expected since the gravitational field will become dom- 
inated by the collective effect of many particles acting on 
any given particle. 



IV. DISCUSSION 

In this article we have studied in detail the early time 
evolution of infinite self-gravitating shuffled lattices, as 
well as the limiting case given by the Poisson distribu- 
tion. We have shown that a very good description of the 
evolution of two point correlations in this phase is given 
by a simple approximation in which the force on particles 
abruptly switches from that given by the PLT approxi- 
mation developed in 0, ^ to the force due only to NN 
particles. Further in the first phase the system evolves at 
small scales to always resemble closely the Poisson dis- 
tribution, explaining the universality of the form of the 
non-linear correlation function which emerges. We have 
noted, however, that this universality will not extend to 
SL initial conditions with arbitrarily small initial shuf- 
fling. In this limit effects come into play in the very long 
first (PLT) phase leading to a strongly correlated evolu- 
tion at all scales, which is different from (and unrelated 
to) that in the Poisson distribution. 

We have thus given, for this specific class of initial 
conditions, an explanation of the non-linear correlations 
which emerge at early times. As underlined in |3|, Q] this 
non-linear correlation function coincides with that which 
is observed at later times, when the system manifests 
a simple spatio-temporal scaling (or "self-similarity"). 
Thus the model appears to explain this asymptotic form 
of the non-linear correlations. As described in [Q this 
can be understood also in the following way: these non- 
linear correlations in the system evolving at any later 
time can be well approximated by those in an evolved 
coarse-grained "daughter" distribution. In the latter the 
system may be in the early time phase studied here, while 
the original distribution is not. The non-linear correla- 
tion functions of the two systems nevertheless coincide. 

Our results are of relevance to simulations of structure 
formation in the universe in cosmology. In this case the 
goal of numerical simulation is to recover the non-linear 
correlations in the Vlasov-Poisson (VP) limit of the evo- 
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lution of a self-gravitating iV-body system. These initial 
conditions are different to those used here — with ini- 
tially correlated displacements and an expanding space 
- but, as shown in ||, ||, the evolution is qualitatively 
similar to our simpler case. If the same kind of model 
can be used to explain non-linear correlations in an early 
time regime in this case then these correlations clearly 
are not described by the VP limit: the forces due to NN 
particles are neglected in this (mean field) limit, while in 
this model they are the dominant ones. In this context 
this would mean that results in this regime would need 
to be discarded as unphysical (i.e. not representing the 
required physical limit, but just a numerical effect aris- 
ing from the method of discretization). Further, if the 
form of the asymptotic non-linear correlation function is 
really determined by this early time evolution, these ef- 
fects of discreteness are then important at all times and 
the system never represents the VP limit. This resem- 
blance of the early time and asymptotic correlation func- 
tion may, however, be a simple coincidence. The evolu- 
tion at longer times may then indeed be representative 
of the VP limit while discreteness effects are important 
only at early times. 

We will study these issues further in future work. 
Firstly the application of our simple two phase descrip- 
tion to cosmological simulations should be investigated. 
We expect the expansion of space to change nothing qual- 
itatively in our model: the description of the PLT phase 
is qualitatively unchanged S, and NN domination dur- 



ing the formation of the first non-linear structures in 
such simulations has been explicitly shown in fusf . The 
fact that the typical initial conditions of such simula- 
tions have much more long wavelength power than those 
considered here may, however, be important [typically 
one has P(k) ~ k m with m < — 1, compared to n = 2 
(SL)and n = (Poisson) considered here]. Secondly, as 
we have discussed at some length in the conclusions of 
jl| , it would be instructive to study numerically the rela- 
tion of the asymptotic regime to the early time regime by 
completely modifying the former using a large smoothing 
in the force, i.e., a smoothing scale e ^ I. In this case the 
dynamics we have described, between NN particles, will 
not occur. However at longer times the system should 
still evolve the same correlations if the VP limit of the 
system is indeed that represented by the simulations with 
e < I. 
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